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Abstract 

The stabilizing effects of local enrichment are revisited. Diffusively coupled host-parasitoid and 
predator-prey metapopulations are shown to admit a stable fixed point, limit cycle or stable torus 
with a rich bifurcation structure. A linear toy model that yields many of the basic qualitative 
features of this system is presented. The further nonlinear complications are analyzed in the 
framework of the marginally stable Lotka-Volterra model, and the continuous time analog of the 
unstable, host-parasitoid Nicholson-Bailey model. The dependence of the results on the migration 
rate and level of spatial variations is examined, and the possibility of "nonlocal" effect of enrich- 
ment, where local enrichment induces stable oscillations at a distance, is studied. A simple method 
for basic estimation of the relative importance of this effect in experimental systems is presented 
and exemplified. 

PACS numbers: 87.17.Aa,05.45.Yv,87.17.Ee,82.40.Np 
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I. INTRODUCTION 



Spatially extended ecological systems, dispersal and metapopulation dynamics have at- 
tracted a lot of interest in recent years |l( ■ m particular, the stability of an ecological system 
that contains many species interacting via competition, predation and symbiosis, and subject 
to environmental and demographic noise, poses an interesting mathematical and ecological 
puzzle. The role of dispersal among sites, migration amongst habitat patches, "rescue" 
of a habitat by dispersal from other locations and so on, have been recognized as crucial 
ingredients in the stabilization mechanism. 

This work concentrated on one of the basic ecological processes - the dynamics of victim- 
exploiter (predator-prey, host-parasitoid) populations. In the "mean field" limit, i.e., where 
the chance for a predation event is fixed for any exploiter- victim pair, the simplest mathemat- 
ical models are those of Lotka and Volterra 0; 3] for predator-prey, and Nicholson-Bailey 
4| for host-parasitoid populations. Both models do not allow for an attractive manifold 
(such as fixed point, limit cycle or strange attractor). In fact, the coexistence point is either 
marginally stable (Lotka- Volterra like) or unstable (with the Nicholson-Bailey host para- 
sitoid system as a prototype). In both cases one may expect to see extinction of (at least) 
one of the species after a short time; if the system is marginally stable the noise will drive it 
to extinction, while for an unstable manifold the system gravitates perpetually towards the 
edge of extinction. 

How do victim-exploiter systems persist in nature for millions of years? One may suggest 
that the simple models are wrong, and assert that any realistic mathematical description of 
a victim-exploiter system should be " dissipative" , supporting an attractive manifold. There 
are many modifications of the basic models that achieve these goals (e.g., by taking into 
account the finite carrying capacity of the environment or adding delayed response), and 
the observed population oscillations may be a result of noise perturbing an attractive fixed 
point [5| . However, since the work of Gause 6J through the classic experiments of Pimentel 
et al. [7], Luckinbill [s] and Huffaker [3], it was known that small sized predator prey (or 
hosts and parasites) systems reach extinction in experimental time scales. The reason for 
this has become clear in the last decade, due to the experiments of Holyoak and Lawler 10) , 



Kerr et al. 



11 



12] and Ellner et al. 131 ] ■ In all of these experiments the same system goes 



extinct rapidly in the well-mixed limit while persisting way above the experimental time (up 
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to h undreds of generations) when the population is spatially segregated. It follows that in 
many cases, and perhaps even generically, victim-exploiter systems are unstable in the well 
mixed limit, and acquire their stability due to its spatial structure. 

What exactly stabilizes the spatial structure? a few answers have already been suggested. 
Some are related to the effect of spatial heterogeneity or spatio-temporal fluctuations, while 
others stress more generic mechanisms that may yield stability even if the space is perfectly 
homogenous. In this paper we focus on the effect of spatial heterogeneity (e.g., different 
growth rates on different spatial patches) as a stabilizer. 

The observation that diffusive coupling between patches may stabilize otherwise unstable 
dynamics, or may result in convergence to a focus instead of a limit cycle, has been made in 



few disciplines independently. In ecology, Murdoch and Oaten 



14j suggested that dispersal 



between Lotka-Voltera patches with spatial variability may stabilize the fixed point. Subse- 
quent works by Crowley [l5|, Ives \l6\. Murdoch et. al. 17] and Taylor lgj], extended this 
basic idea to include multi-patch systems, the effects of parasitoid aggregation, differences in 
diffusion parameters, density different migration and other complications that may occur in 
realistic systems. In chemistry, this stabilization is known alternatively as "oscillator death" 
and was observed by Bar-Eli [19| in the context of coupled chemical oscillators. That basic 
idea has since been applied to other diffusively coupled chemical systems, such as neural 
[20I and calcium oscillations [21I . Mathematically speaking, the stabilizing effect of diffusive 
coupling between two sites on a single species, extinction prone chaotic system has been 
considered by Gyllenberg et. al. [22]. These authors pointed out the "salvage effect", where 
the existence of a sink may stabilize the population on the source habitat; this effect also 
appears in the systems considered below. 

Although the basic phenomenon is known, it turns out that the system of diffusively 
coupled unstable oscillations is quite rich and may reveal many interesting features beyond 
the stabilization of a fixed point. In this paper we analyze a very simple case of a single 
enriched site (where the prey, or host population flourish) surrounded by a less productive 
environment. We will examine in detail the conditions for the stability of the fixed point 
and the asymptotic behavior of the Lyapunov exponent, consider the possibility for the 
appearance of a limit cycle, and discuss the associated bifurcations. The most interesting 
phenomenon observed relates to the spatial population profile: here the effect of localized 
enrichment may yield either local or nonlocal changes, including the emergence of oscillations 
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far from the location of enrichment, incommensurate oscillations, etc. 

This paper is organized as follows: In the next section a linear toy model is presented and a 
few basic insights are derived. The effect of nonlinearities is emphasized for a predator-prey, 
marginally stable model in the third section. Although marginal stability is not a robust 
feature of a dynamical system and should not be considered seriously as the underlying 
dynamic of ecological processes, its consideration enables the attainment of basic knowledge 
regarding the effect of nonlinearity on stability, as will explained later on. The fourth section 
analyzes a Nicholson-Bailey like (unstable) system and discusses interesting nonlocal effects 
of enrichment. Finally, our results are discussed in view of the "paradox of enrichment" 
and possible experimental tests are suggested. 



II. A TOY MODEL: COUPLED LINEAR OSCILLATORS 

Let us present a linear model in order to exemplify some features of the fixed point sta- 
bility (local analysis) for non-identical victim-exploiter patches. We consider two diffusively 
coupled harmonic oscillators (with the possibility of repulsive term with a positive constant 
a): 



— = u 1 y 1 + D{x 2 -ari) + axi 
dx 2 

oj 2 y 2 + D(xi - x 2 ) + ax 2 



dt 

= -u 1 x 1 + D(y 2 - y x ) + ay x (1) 



dyi 



dt 

— = -uj 2 x 2 + D(y 1 - y 2 ) + ay 2 . 

As this system is linear, it may be diagonalized around the (only) fixed point at zero. 
When a = 0, the Lyapunov exponent T for that fixed point turns out to be negative as long 
as \8\ = uj 2 — lo\ 7^ for any D, and approaches zero (marginal stability) if the dispersal is 
very small (no connection between oscillators), very large (single oscillator limit) or if the 
system is homogenous (8 — > 0). The Lyapunov exponent is given by: 



T = a + Re 



-D + 1/2 \l AD 2 - AujxS + 2 J- (5 + 2wi) 2 (-5 2 + AD 2 ) - 4cji 2 - 25 2 



(2) 
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and its typical behavior is illustrated in Figure [TJ Linearity implies that if the fixed point 
is stable it is also globally attractive. The parametric dependence is characterized by the 
following properties: 

• Without loss of generality lo\ may be scaled to unity by rescaling of the time. Thus 
the stability is determined by three parameters: migration rate, repulsion (a) and the 
desynchronization term 5. 

• T is a nonmonotonic function of the migration rate; close to zero migration, (T — a) 
vanishes linearly with D, while for large diffusion it decays like 1/D. The optimal 
dispersal (given that other parameters hold fixed) is D = \5\/2; in which case T = 



• The only effect of a is a rigid upward shift of T, as demonstrated in Figure HJ where 
the dashed line indicates the border between the stable and the unstable regime. 

• For fixed migration an increase of 5 always helps to stabilize the system, but the effect 
saturates at 5 = 2D. Accordingly, for any a, there is a critical diffusion below which 
the system turns to be unstable, independent of the level of heterogeneity. 

The limitations of this toy model are related to its linearity. In particular, once the 
fixed point at zero is unstable, the oscillation amplitude will grow unboundedly, and no 
other attractive manifold is allowed. This should not be the case if the coupled oscillators 
are nonlinear. One possible mechanism that may be responsible for an attractive manifold 
becomes clear if one looks at the system ([1]) in polar coordinates, where = xf + yf, 
9i = arctan(yi/xi) (i — 1, 2). In this representation the overall phase 6i + 62 decouples and 
the phase space turns out to be three dimensional, 



Where the phase desynchronization is = 62 — 61, the amplitude desynchronization is 
r = r 2 — ri, and the homogenous manifold is one dimensional R = r% + r 2 . As implied from 
([[]), the dynamic is either a flow towards an attractive fixed point at R = 0, or an unbounded 



a- \8\/2. 




(3) 
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FIG. 1: The Lyapunov exponent of the system (TT|), plotted vs. the migration rate D. Parameters 
are ui = 1. and u>2 = 1.2, 1.4, 1.6, 1.8, as indicated by the legend; a = 0.2. Note that the system 
attains its maximal stability as D = 5/2, as predicted. The inset shows the same lines for a = 
in a log-log plot, to emphasize the asymptotic behavior D and 1/D. 

growth of R, depending on the relation between the stabilizing desynchronization factor 5 
and the repulsion a. If the system is nonlinear, however, the angular velocity u is generally 



r dependent 



241 ] . This implies that, although 5 is too small close to zero, the "effective 5" 



may be larger far from the fixed point, leading to desynchronization and stabilization. In 
such a case one may expect a limit cycle in the nonlinear system. This, in fact, actually 
occurs, as will be demonstrated in the next sections. 

Another limitation of the linear model is the fact that the location of the fixed point 
(at zero) is independent of the migration rate. Equations describing population dynamics 
support a nontrivial coexistence fixed point as well, and diffusion may alter not only its 
stability but also its location. In the next section we will show that the results for the 
asymptotic behavior of the Lyapunov exponent, as well as the saturation of T for large S, 
may not hold in more complicated nonlinear dynamics. 
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III. MARGINALLY STABLE LOTKA-VOLTERRA DYNAMICS 



Let us now take one step towards more realistic victim-exploiter systems, and consider the 
Lotka-Volterra continuous time dynamics, denoting the density of the predator population 
by a, and the prey density by b. On a single patch (say, when the population is well 
mixed, and the probability of encounter and predation is equal for any two individuals in 
the population) the pure LV system is marginally stable; there is a coexistence fixed point, 
but the system may also oscillate with any amplitude around this point, and the dynamic is 
determined by the initial conditions. Any type of stochasticity (demographic, environmental, 
random migration) will lead to random wandering of the system between all possible orbits. 
The amplitude of oscillations thus performs some sort of random walk, with the average 
amplitude growing in time until extinction 24]. 



To consider the stabilizing effect of spatial differences in the environment, we will focus 
our attention on two typical situations: the case of two coupled patches, where the basic 
stability properties are to be demonstrated, and the case of a one dimensional " chain" (with 
N spatial patches and periodic boundary conditions) where the spatial population profile is 
examined. 

We assume that predator-prey dynamics have the same parameters on all sites, the only 
exception is the "zero" site where a, the growth rate of the prey, is larger, i.e.,: 

J a ° n = t*\ 

°n = < (4) 

I 0i else. 

The corresponding Lotka-Volterra equations for a one dimensional array are: 

-fMa n + \Aa n b n + D A (—2a n + a n+ i + a n _i) 



dt 

= a n b n - \ B a n b n + D B (-2b n + b n+1 + 6 n _i), (5) 

where /i is the predator death rate, A is the predation rate and Da and D B are the hopping 
rates of animals from one spatial patch to the other. Note that one can use the rescaling of 
the densities in order to take \a = A# = 1 and set // = 1 using the rescaling of time; this is 
the parametrization used hereon. Moreover, we focus here on the case of equal diffusivities 
D A = D B = D. 

As our system consists of one special site connected to a "reservoir", it turns out that 
the character of the reservoir is also of importance. A distinction should be made between 



the case of a single "oasis" coupled to a desert, (i.e., where o\ < so the linear growth rate 
of both predator and prey is negative on all the patches except one) and the case where 
Co > 0"i > 0, where all sites are "active". Hereon, the first case will be regarded as an 
"oasis-desert" situation (OD), where the case o~i > is named the poor-rich (PR) scenario. 

To begin, let us consider the two-patch case and examine the stability properties of the 
system. As the LV is marginally stable, one may expect that any spatial heterogeneity will 
yield stability. This is, indeed, the situation, as demonstrated in Figures [2j The general 
structure is very close to the linear model predictions with a few exceptions. First, the D 
dependence of the Lyapunov exponent at small dispersal holds only in the PR case, while 
the OD scenario is characterized by D 2 asymptotic. Second, for the large D asymptotic, 
T approaches zero like 1/D 3 instead of 1/D in the linear model. Both asymptotes may be 
obtained analytically as explained in appendix A. Note also the rightward migration of the 
optimal diffusion as the difference cr 2 — o\ increases; this is the analog of the 8/2 dependence 
of the linear model. As previously explained, the effects of nonlinearity have to do with 
the fact that the location of the fixed point itself depends on the migration rate. Beyond 
the two-patch limit, the stability properties are qualitatively the same, as demonstrated in 
Figure [31 

For the extended system one may consider not only the stability properties but also the 
spatial profile of the population densities. Let us consider first the oasis-desert scenario. 
Numerically solving the spatial configuration for (jSJ) with the spatial heterogeneity defined 
by (HI), one gets the colony profile presented in Figure HI Clearly, deep in the desert the 
populations of both prey and predator are small, so the nonlinear (predation) terms in (J5J) are 
negligible. Accordingly, each of the species follows asymptotically the profile of a logistically 
growing population around an oasis [l|; Q; 26], and the density decays like exp(— a\\x\/D) 
for the prey, and like exp(— n\x\/D) for the predator. If we adopt the definition of the size 
of a colony as the length scale for which the population density is equal to some constant 
(this reflects the threshold introduced by the discreteness of the individuals), this size (up 
to logarithmic corrections) is 1§Db/(J\ for the prey and IqDa/ H for the predator, where lo is 
the lattice spacing. 

This result has two implications. First, in case of multiple oases, the interesting quantity 
is the typical distance between them in units of the colony size; if this is a large number, 
one should consider a system of two different patches, while a small number indicates strong 
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FIG. 2: r vs. D for the two patch system. The Lyapunov exponent is calculated by numerical 
diagonalization of the stability matrix obtained by linearizing Eqs. © around the fixed point, 
where the prey fecundity is do = 0.05 and o\ = 0.02, 0, and —0.01, as indicated by the legends. 
In all cases, the coexistence fixed point is stable for any finite migration. The log-log plot (inset) 
emphasizes the asymptotic behavior in the limits of small and large dispersal: if a\ > 0, the linear 
model predictions still hold and T ~ D at small D, while for a\ < (coupling of an oasis to a 
desert) F ~ D 2 . For large migrations, on the other hand, in all cases the asymptote differs from 
the linear model predictions, and T ~ 1/D 3 . 

mixing. The chance for a "rescue effect" (where, due to noise, the population on one patch 
goes extinct until the "rescue" by a rare immigrant from another habitat) may be very 
different for the victim and for the exploiter, depending on their death- diffusion ratio in 
the desert. Second, as demonstrated in Figure HI the spatial decay of the colony may give 
false hints regarding the density at the oasis: in the example presented here, the predator 
population is rarely far away, but the predator abundance is much larger than that of the 
prey on the oasis. While the linear death rate of the prey is constant along the system and 
determines the exponential decay, its population on the oasis is dictated by the prey growth 
rate ctq and is much higher. 
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FIG. 3: r vs. ln{D) for two, four, six and eight patch system. The Lyapunov exponent is the result 
of a numerical diagonalization of the stability matrix obtained by linearizing Eqs. ([5]) around the 
fixed point. 

The effect of local enrichment becomes even more interesting in the poor-rich case, as one 
sees in Figure [51 An enrichment of the prey growth rate (say, by increasing the food supply) 
on a site leads to depletion of the prey in its neighborhood, while leaving the population on 
the rich site above its pre enrichment level. Surprisingly, the behavior of the prey density is 
paradoxical, but the "paradox" is nonlocal. 

IV. UNSTABLE SYSTEMS: NICHOLSON-BAILEY LIKE MODEL 

The original mathematical description of a host-parasitoid system was formulated by 
Nicholson and Bailey [4j. If the host density is given by H and the parasitoid is P, Nicholson 
and Bailey map for nonoverlapping generations is: 



qH t e 



-zPt 




(6) 
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FIG. 4: Colony profile for a diffusively coupled oasis-desert system (ctq = 2.5, o\ = —0.25, D = 1) 
Note that the ratio between densities at the oasis differs from that ratio at the desert. 



where q is the fecundity of the host in the absence of the parasitoid, z is the infectivity (the 
chance of a single parasitoid to infect the host) and c is the number of parasites produced 
by a single infection. 

The essential feature of the NB model is its instability - it allows ever growing oscilla- 
tions around its coexistence fixed point. However, the model does not allow explicitly for 
extinction, as P and H are positive definite along the process. As the oscillations grow, the 
minimal distance from the extinction phase for each of the species decays rapidly. Hence, 
for any realistic system, where small noise or the discreteness of individuals are taken into 
account, the system flows deterministically to the extinction phase. 

One technical characteristic of the NB dynamics makes its analysis problematic in the 
context of the current discussion. Eqs. ([6]) assume nonoverlapping generations, and the cor- 
responding map, when extended to diffusively coupled spatial domains, supports attractive 
periodic orbits with finite basin of attraction 22J. Even in the absence of spatial inho- 
mogeneities some initial conditions are trapped in the periodic orbits, and it is hard to 
distinguish between this effect and stabilization due to spatial differences. In order to focus 
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FIG. 5: Colony profile for a diffusively coupled poor-rich system (<ro = 43.5, o\ = 27, ,D = 
0.5). Without the diffusive coupling, each site supports a (marginally stable) fixed point for the 
predator (dashed black) and for the prey (dashed red). In the diffusively induced stable fixed point, 
the predator migration increases the density of predators in the poor sites close to the rich one, 
thus depleting the local prey population, and the overall prey profile admits a minimum in the 
neighborhood of the rich site. 

the discussion on the consequences of local enrichment, we will switch here to continuous 
time dynamics that imitates the relevant features of the NB model. For the sake of nonlinear 
dynamics analysis considered here, our model is equivalent to the small z limit of the NB 
system, where the map may be approximated by continuous time equations, as emphasized 
in figure El 

Let us assume, thus, that in a predator-prey system the encounter between a single 
predator and prey may result in predation, but the predation probability increases if two 
predators encounter a single prey (i.e., there is a predation "Allee effect"). We "enrich" the 
LV system with additional terms that correspond to this ally predation process: 
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FIG. 6: Phase portrait, calculated by numerical integration of Eqs. shows the main character- 
istics of the NB continuous analogue: the instability of the fixed point, the exponential growth of 
the oscillations, and the fact that the trajectory cannot cross the boundaries of zero population 
for either predator or prey. The corresponding trajectories of the map ([6]) for host-parasitoid are 
shown in the inset. 



da , 2l 

— = — ua + Aaclo + aa o 
at 

— = ob — Xsab — aa 2 b. (7) 
dt w 

The resulting equations admit a coexistence fixed point at a = (A — A)/2a, b = 2/i/(A + 
A), where A = VA 2 + 4aa. This fixed point is unstable, with a positive local Lyapunov 
exponent. The amplitude of oscillation increases rapidly with time, but the system cannot 
reach the extinction phase unless some noise is presented. 

To include spatial structure, we switch to a nondimensionalized form of the diffusively 
coupled patch system, where the definition of constants is identical with Eq. [5j 
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da. 



if 



— a n + a n b n + aa: 



+ D A (-2a n + a n+1 + a n _i) 



dt 



0" n &„ - a n b n - aa^6 n + D B (-2b n + b n+1 + fe„_i). 



(8) 



Again, one should expect to recover the zero dimensional dynamics as D approaches zero 
and infinity. In between, the resulting dynamics turn out to be quite rich. As in the LV 
case, we will check the two-patch system to get the general picture for the stability of the 
system, then look at a multi-patch example to see the spatial profile. 

The linear model admits two options: in the presence of spacial variation the unstable 
system may remain unstable (in which case it flows to extinction exponentially fast) or 
admit (for intermediate diffusivities) an attractive fixed point. As exemplified in Figured 
the nonlinear model (jSJ) turns out to admit another option: while the fixed point may be 
unstable, the system does not flow to extinction but to an attractive manifold like a limit 
cycle. This gives rise to many complications in the system behavior, including various types 
of bifurcations, bistability and hysteresis. Partial analysis of the bifurcations and the phase 
portraits for the unstable case is presented in Appendix B. 

As in the marginally stable LV case, one observes a nonlocal effect of enrichment on the 
spatial profile. In Figure [8] a few possibilities are demonstrated. It may happen that nutrient 
enrichment at a single site causes the other parts of the system to change periodically, while 
the rich site population is almost fixed. Increasing the spatial variation, the amplitude of 
oscillations on the rich site grow (Figure EJ), and at the same time the phase locking with 
the "wings" is lost. At the end, the oscillations at the center and on the wings become 
incommensurate, giving rise to an attractive torus, as demonstrated in Figure OH In the 
oasis-desert case the effect is less dramatic, and the oscillations appear first in the enriched 
region (Figure ITT1). 

V. DISCUSSION 

The effect of quenched spatial variations, where the habitat is made of local patches 
connected by dispersal, is stabilizing. Besides the survey of the possible states of this 
dynamical system and the asymptotical behavior for large and small migration rates, the 
nonlocal effect of enrichment emerges from the numerics as a new feature. In particular, 
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FIG. 7: Stability diagram for the coexistence fixed point of the two-patch Nicholson-Bailey like 
system, Eqs. (jSJ). Parameters are o"o = 1, a\ = —0.1, and a = 0.1 (low) a = 0.4 (middle) and 
a = 1 (upper line). Only the low a allows for a finite region of dispersal rate in which a stable 
coexistence fixed point is supported. However, even in the unstable region the system may flow 
into a limit cycle or other attractive manifolds, see Appendix B. The inset shows a stable orbit 
(projected on the homogenous manifold) for the system parameters at the point indicated by an 
arrow. 



it turns out that local enrichment may stabilize otherwise unstable system, and that the 
limit cycle may involve periods of large amplitude far from the enriched site, while the rich 
site itself involves cycles of negligible amplitude. Increasing the spatial heterogeneity, one 
observes an increase in oscillation amplitude on the rich site, while its neighborhood remains 
relatively calm. The oscillations at the enriched site may be incommensurate with the period 
of oscillations at the "wings", such that the system flows into a stable torus instead of a 
limit cycle. 



This observation may solve the "paradox of enrichment" 23j in some cases. Many victim- 
exploiter models predict an increase in oscillation amplitude as a result of increased prey 
growth rate or carrying capacity. The empirical support of this prediction, however, is quite 
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FIG. 8: Population profile for a 16 site system with periodic boundary conditions (a few snapshots 
taken along the evolution of the system) for the prey (red) and the predator (black). The two upper 
panels correspond to a = 0.02, an = 5, a\ = 2.25 and D = 0.5. One notices that the enrichment 
almost fixes the population on the rich (zero) site at the middle, while in the far zone (the "wings") 
the system oscillates. The lower panels show what happens as the spatial variation is increased: 
here un = 21.5 and all other parameters are the same. Due to the increase of heterogeneity, cycles 
also developed on the rich site, while in the sites close to the origin, their amplitude remains 
negligible. 



29; 



30] . From the results here 



limited [28( and in many experiments the effect is absent 
one concludes that in some cases, small spatial fluctuations in the enrichment, e.g., small 
inhomogeneities in the concentration of food, sunlight or other resources, may stabilize the 
system and avoid the "paradoxical" behavior. 

What degree of enrichment uniformity is required in order to observe the paradoxical 
behavior? Such a question may be analyzed very easily using the toy model presented 
here as a basic tool for parameter estimation. As stressed above, only three parameters 
control the system close to the fixed point: the dispersal rate D, the repulsion a and the 
desynchronization factor 5. The threshold between a " uniform like" and heterogenous system 
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FIG. 9: Prey (red) and predator (black) populations vs. time for crn = 5 at the "wings" (far 
zone - see figure EJ (panel la) and at the rich point (panel lb; the prey at the enriched point is 
not presented as its oscillations are relatively small). Panel 2a shows the prey and predator at 
the wings for an = 21.5, while in 2b the oscillations at the rich point are graphed for the same 
heterogeneity. 

appears at D ~ 5, i.e., when the migration rate in/out of a patch (here a patch is defined as 
a region where the resource distribution is uniform) is close to the desynchronazation rate, 
given by the difference in oscillation period between patches. If diffusion is much larger, the 
system should be considered as uniform. Much smaller diffusion corresponds to the situation 
of almost independent patches. 

For an experimental system, one should try to gain a rough estimate of the above parame- 
ters in order to ascertain the level of enrichment needed to yield a strong enough effect. Here 
we exemplify these considerations using two experiments: one on a predator-prey system 
and the other host-parasitoid. 



In the experiment of Kerr et al 12], E-coli bacteria is the host and its viral pathogen, 
T4 coliphage, is the parasite. In the well-mixed case, the system survives about 24h before 
the bacteria undergoes extinction; this corresponds to a ~ 10~ 5 seconds. The migration in 
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Enriched "Wings" "Phase" 




FIG. 10: Phase portrait of the populations (predator vs. prey) for the cases considered in figures 
[8] and EJ The do = 5 case yields the upper panels, where la is the phase portrait at a point on 
the wing, lb is on the enriched site (not the low amplitude of oscillations) and in lc the prey 
at the center is graphed vs. the predator on the wings. 2a and 2b are the same graphs for the 
an = 21.5 case, and the fact that the oscillations are incommensurate manifests itself in panel 2c. 
In all subfigures the y axis is the predator density and the x axis is the prey density. 

that experiment is controlled by the experimentalists, as a robot taking biological material 
between otherwise disconnected patches. The authors present simulations, based on empir- 
ically calibrated cellular automata, that predict oscillations with a time scale of about 10 
days, so uj\ is about 10~ 5 , i.e., close to a. If one chooses a day as the unit time, uj ~ a ~ 1. 
From the linear system analysis, it seems that this system should be coupled to another one 
with co>2 ~ 2co>i in order to observe stabilization due to spatial heterogeneity. 

Another example is the experiment of Holyoak and Lawler [lol |. where the predaceous 
ciliate Didinium nasutom is feeding on the bacterivorous ciliate Colpidium cf. Striatum. 
The diffusion constant for these small (length of order 0.1-0.05 mm) creatures depends on 
the level of water turbulence and on the size of the tubes connecting different microcosms. 
The extinction time in the well mixed case is about 70 days (a ~ 5 * 10~ 6 seconds) and U\ 
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FIG. 11: Prey (red) and predator (black) populations vs. time for a® = 21.5, o\ = —0.5, a = 0.02 
and D = 0.5. Here the oasis (at zero) is coupled to a desert, and hence the oscillations take place 
only at the enriched location. 

is about 10 days. In the appropriate dispersal, one may observe stabilization if the system 
is coupled to another set of microcosms with 5 = 002 — u)\ ~ 0.05, i.e., relatively small 
enrichment will allow for stabilization. 

We acknowledge helpful discussions with Marcel Holyoak. This work was supported by 
the EU 6th framework C03 pathfinder and DAPHNet. 

VI. APPENDIX A 

In this appendix we discuss the dependence of the system's stability (its Lyapunov expo- 
nent) on the rate of animal dispersal. As explained above, one should recover the mean field 
results in the case of no migration, D = 0, and as the migration rate approaches infinity. 
For the Lotka-Voltera, marginally stable case, the Lyapunov exponent vanishes in these two 
limits, while in between it is finite and negative. We now want to explore the two limits 
more carefully and attain the asymptotic dependence of the Lyapunov exponent in terms of 
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the system parameters. 

To begin, let us consider the simplest spatially explicit CclSG, 1.6., db two-patch LV system. 
The dynamics in such a case are four dimensional and described by: 

— = -ai + ai&i + D(a 2 - aj 

— = — a 2 + a 2 b 2 + D{ai - a 2 ) 

— 1 = a x bi - aifei + D(b 2 - 6 X ) 
db 

-l = C r 2 b 2 -a 2 b 2 + D(b 1 -b 2 ), (9) 

where the growth parameter a is heterogenous. Clearly, as D approaches infinity the differ- 
ences between patches, p = a\ — a 2 and 6 = b\ — b 2 approach zero, so it is useful to rotate the 
coordinate system in order to separate between the homogenous manifold (A = (a\ + a 2 )/2 
and B — (b\ + b 2 )/2) and the p — 9 manifold. Defining, now, e = l/D, and assuming that 
A = Aq + eA\... (and the same for B) while p = tp\ + t 2 p 2 ... (and the same for 6), one can 
solve for the fixed point of the set of equations order by order in e up to, say, order n. 
This solution is then plugged into the stability matrix, and the roots of the characteristic 
polynomial may then be found (again, order by order in e) up to 0(n). With that, one 
finds that the leading contribution to the Lyapunov exponent (the first order for which the 
root of the characteristic polynomial admits a real part) is n = 3. The asymptotics of the 
Lyapunov exponent for large D turns out to be: 

(ai - a 2 ) 2 ((Ji + 02) _ 



64L> 3 

Close to D — 0, on the other hand, one should use another technique. Solving 01,02,61 
and b 2 to the n-th order in D, these have to be plugged into the stability matrix, from which 
the characteristic polynomial is extracted. For D = 0, two (in the oasis-desert case) or all 
four (in the poor-rich case) of the eigenvalues are purely imaginary. The corrections to these 
imaginary eigenvalues are written as a power series in D, and one finds that the real part 
of the slowly decaying eigenfunction (i.e., the Lyapunov exponent) is O(D) in the poor-rich 
case, but only 0(D 2 ) for the oasis-desert situation. For a 2 < (oasis desert situation) one 
gets: 

r _ Q-2 3 + 0-2 0-1 + O-l 2 + 0-1 D 2 , . 

O-2 (CXI 0- 2 2 + 0- 2 2 + ai 2 + O-i) 2 ' [ } 
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While in the poor-rich case for o\ > u% the exponent is: 

T = -^±^D. (12) 
2o"i 

This result does not hold when the system approaches the uniform limit o\ = 02, as it has 
been obtained perturbatively, neglecting the degeneracy appearing for these parameters. 



VII. APPENDIX B 



In this appendix we consider the various equilibrium states of the unstable dynamics 
(J7J) on a heterogenous environment (jSJ). The 4d system is quite complicated and allows for 
many types of bifurcations. In the following, a general sketch of the main types of dynamical 
behavior is given. 

Let us take a look at Figure [121 which describes the Lyapunov exponent of the coexistence 
fixed point vs. D. In the region of intermediate migration (between 2 and 3), the coexistence 
point becomes attractive due to the heterogeneity, similar to the LV system. However, at 2 
the system undergoes a supercritical Hopf bifurcation, and a globally attractive limit cycle 
appears. As the diffusion constant approaches zero the radius of this limit cycle diverges, 
and finally it losses its stability via homoclinic bifurcation, and the cycle period diverges 
logarithmically as shown in the inset to Figure [T2l 

There exists, however, another set of parameters, where the coexistence fixed point does 
not acquire stability for any diffusion rate, as demonstrated in Figure [T3l It turns out that, 
in this case, the limit cycle born via homoclinic bifurcation at D — loses its stability 
through an infinite period bifurcation at larger dispersal values, leaving an unstable system 
with oscillations that grow unboundedly until extinction. The stable cycle reemerges at even 
higher values of D, while now its size grows as diffusion increases. It should be noted that 
in some cases (like the parameter range that corresponds to the middle curve in Figure [7]) 
we have observed a limit cycle behavior for any value of D. 
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FIG. 12: Nicholson Bailey system ([8]), where the parameters (specified in Figure [7]) allow for a 
locally stable coexistence point. The bifurcation takes place at the points indicated by 2 and 3, 
and between 2 and 1 a limit cycle exists. This limit cycle annihilates with the saddle via homoclinic 
bifurcation, characterized by a logarithmically diverging period of oscillations, as indicated in the 
inset. For D values right above the upper point 3 we have observed trajectories that converge to 
a limit cycle, but this cycle disappears close to that point and leaves an unstable system. 
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FIG. 13: Nicholson-Bailey like system with large a. Here there is no stable fixed point for any value 
of migration, yet it may support a limit cycle. This limit cycle exists only to the left of the dot 
indicated by the arrow. At this point the limit cycle disappears in an infinite period bifurcation. 
The time spent by the system in its unstable phase (to the right of the indicated point) close to 
"ghost" of the limit cycle diverges close to the bifurcation like \j\fe, where e = D — D c is the 
distance from the bifurcation point (inset). 
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